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1. Introduction 

Our understanding of complex physical systems is often aided by the construction of 
simple mathematical models. A prime example from glass theory is the trap model 
of Bouchaud |1|, 0, which describes a system's relaxation purely in terms of activated 
processes between potential energy wells of various depths. With only a few assumptions 
the model manages to reproduce many of the phenomena commonly associated with 
glasses, such as a crossover to non-equilibrium behaviour as the temperature is lowered 
and aging in the low-temperature regime. Furthermore the model is sufficiently generic 
to apply to structural as well as spin glasses, and has recently been extended to the 
study of rheology in glassy systems 

However, there are a number of spin-glass models that exhibit a continuous phase 
transition, or critical point, in the limit of zero temperature, quite unlike the current 
realisation of Bouchaud's trap model. These include Stein and Newman's model, which 
maps onto invasion percolation at a temperature T = the two dimensional Ising 

spin glass P]; and the Sherrington-Kirkpatrick model under infinitesimal variations of 



the external magnetic field [10 



In this paper we introduce a modified version of Bouchaud's trap model which 
incorporates a form of interaction between different subsystems. It is shown that in the 
limit of zero temperature this model maps onto a self-organised critical system known 



as the Bak-Sneppen model [|Tl|, |T^, so T = 0^ is a critical point in all dimensions, 
including the infinite-dimensional mean field model. Although the Bak-Sneppen model 
was originally devised to explain the pattern of evolutionary bursts in the fossil record 
known as 'punctuated equilibria,' it was justified in terms of activated processes over 
(fitness) barriers and so its biological interpretation is purely semantical. 

This paper is arranged as follows. In section ^, Bouchaud's model is briefiy reviewed 
before the new model is defined in its mean field form. The corresponding master 
equation is solved in the steady state. The relationship with the Bak-Sneppen model 
implies that the model can also be defined on a lattice, and this modification is studied 
numerically in section ^. The simulations were also employed to check the analytical 
predictions of the mean field equations. Finally, our work is summarised in section 0. 



2. Definition of the mean field model and the stationary solution 

In Bouchaud's formalism, the state of the system is represented by a particle (or many 
particles, if the system is treated as a collection of subsystems) moving in configuration 
space over an energy landscape, which is 'rugged' in the sense that most of the particle's 
time is spent trapped in potential wells or 'traps' [|I], 0. Motion within traps is 
not explicitly incorporated; instead, it is simply assumed that whatever dynamics is 
present gives rise an activation probability that depends on the barrier height b and the 
temperature T according to the Arrhenius form u;oe~^^'^, where loq fixes the timescale. 

When a particle is activated it 'hops' to a new trap, where it remains until it is 
again activated by thermal fiuctuations. Thus the state of the system can be described 
by the probability P{b,t) that it is in a trap of depth b at time t (in the many-particle 
interpretation, P{b,t) is the distribution of barriers in all the subsystems). Given that 
the barrier heights are distributed according to some pnor distribution p{b), then P{b, t) 
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evolves according to the master equation 
1 dP{b,t) 



- e-"^ P{h,t) + uj{t) p{h) , (1) 
where uiouiit) is the total rate of hopping at time t, 

/•oo 

uou{t)=UQ e-''/^P(6,t)db. (2) 
Jo 

The first and second terms on the right hand side of (|I]) correspond to hopping out of 
and into traps of depth 6, respectively. For simplicity is has been assumed that p{h) is 
independent of t and T, and barrier heights before and after a hop are uncorrelated. 
When a steady state solution of (|I|) exists, it takes the form 

Poo(&) = lim P(6, t) = e"/^ p{h) , (3) 



using cUoo = lim(_>oo i^{t). For a prior distribution with an exponential tail p{b) ~ e~^/''° 
(which may be generic to glassy systems 0), Poo(^) can only be normalised for T > , 
when it takes the form of a Boltzmann distribution. No normalisable steady state 
solution exists for T < Bq , when the system perpetually evolves into deeper and deeper 
traps. In this manner the model exhibits a phase transition from equilibrium to non- 
equilibrium behaviour at T = 6o , i^. a glass transition. 

2.1. Master equation for strong interactions 

Although the basic trap model described above exhibits a simple and elegant glass 
transition, it has the undesirable feature that the mean barrier height diverges with time 
when T < bo ■ This suggests that there may be some mechanism currently lacking from 
the model which eventually halts its progression into ever deeper traps. For instance, 
it is known that the system will equilibrate at very long times if it can only sample a 
finite number of different configurations [|I|]. Here we argue that introducing a suitable 
form of interaction between the subsystems can achieve a similar effect. 

Bouchaud et al. have considered the effect of allowing each hop to slightly alter the 
barriers of other subsystems , so we are now working strictly in the many-particle 



interpretation. In their scheme, the barrier distribution P{b,t) is perturbed at a rate 
proportional to the overall hopping rate u (t) , giving rise to an extra term in the master 
equation (|l|) of the form 

/■oo 

cuit) / Tib' - b) {P{b\ t)p{b) - P{b, t)p{b')) db' . (4) 
Jo 

For weak interactions, T{b) is narrow and (^) reduces to a diffusion-like term which was 



shown to not alter the essential nature of the glass transition |T3 

Here we consider the opposite extreme of a broad T{b) = z, where z is a constant. 
This corresponds to a form of strong interaction in which every hopping subsystem 
causes z other subsystems to have their barriers assigned entirely new values. This 
erases the memory of the system at a rate zu){t) and, as shall be demonstrated below, 
ultimately removes the glass transition. Subsystems with very large barriers that would 
not normally hop until very late times will now interact instead, which on average lowers 
their barriers towards the mean of p{b) and allows them to hop much earlier than they 
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would in the absence of strong interactions. Thus the interactions introduce a form of 
dynamically generated cut-off into the system. 

It is straightforward to integrate (^) with T(6) = zto give the new master equation, 



1 dP{h,t) 
ujq dt 



-^/^ P(6, t) + uit) p{h) + zuj{t) [pih) - P{b, t)) , (5) 



with the definition of uif) unchanged from (0). The steady state solution of (^) is 
non-Boltzmann, 




Poo(6) = ^ 1 + p{b), (6) 



so Poo{b) ~ p{b) as 6 — > oo for all T, ie. Poo{b) is always normalisable when z > 0. 
This suggests that there is no true finite temperature glass transition, although it has 
yet to be shown that the system always approaches this stationary solution. Numerical 
evidence that this is indeed the case is presented in Section ^ where it is also argued 
that glassy behaviour may persist over intermediate timeframes at low temperatures. 
For the remainder of this section we restrict our attention to analysis of the stationary 
solution (^. 

2.2. Steady state solution 

For the case of a simple exponential prior p(b) = , it is possible to derive an 

exact expression for uJoo for z > 0. Details of this calculation are given in [Appendix A 



here we just quote the final result. In terms of the reduced temperature x = T/bo , the 
expression for is 



sinvrx f 1 ^(--zcUoo)^! „ 

It may at first appear that the factor of sin vrx on the left hand side of (0) would imply 
that Uoo —>■ for integer values of x. However, there is also a simple pole in the sum for 
such X, and in fact these two effects cancel. The equivalent expression to (J^) for integer 
X is 

(ipsEl , ,„ (l + ^) + g '-^^ . (8) 

X{Z+1) \ ZUJoo) ^ x-k 

The numerical plot of (0) and (Q) given in Fig. |l| confirms the monotonic dependence of 
ijJoQ on temperature. 

One might expect that Poo{b) for small z could be expressed as just an 0{z) 
perturbation around the z = solution given in (^. In fact this is only true for 
temperatures x > 2, as we now demonstrate. As z — 0, (|^ can be expanded in powers 
of z and the higher order terms dropped. The identification of the leading order term 
depends upon the value of x. For x > 1, the leading order term is 0{z Uoo) and 

c^oo = (l - ^) ( 1 + 0{z) ) , (9) 
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Figure 1. Plot of looo against the reduced temperature x = T/bo , found from and 
(H) by the method of interval bisection. From top the bottom, the thin lines refer to 
z = 10^^, 10^^, lO^'^, 10"'* and 10"^ respectively. For comparison the z = solution 
uj oc = 1 — 1/a; is plotted as a thick line. 



which is just an 0{z) perturbation around the original solution (|^). This is not true for 
X < 1, when the [zuJooY term is now leading order and 

/ z \ /sin7ra;\ - 

"-"(ttt) {—) ■ (^°' 

This holds as long as the 0{zuJooY term is much larger than the 0{zuJoo) term, which 
corresponds to z^^^ ^ 1. This should be contrasted to the original model, where cUqo is 
not even defined for x < 1. Repeating this procedure for integer x produces the same 
expression as @) for x > 2 and Uoo ~ —(In 2;)"^ for x = 1. 

The consequence of having 2; > becomes apparent at higher temperatures if one 
considers the distribution of barriers as a whole, rather than just uj 00 ■ By substituting 
(§) into the expression for Poo{b) given in (j^) and expanding in powers of z, one readily 
sees that Poo{b) can only be expressed as an 0{z) perturbation around the original 
solution when x > 2; for x < 2, the assumption of a linear expansion breaks down. This 
higher-temperature divergence can be explained by observing that, when z = 0, the 
mean time between hops diverges in the stationary state for x < 2 [H. By contrast, 
when z > the mean time between interactions approaches (zuJoq)^^, which is always 
finite. Thus the interactions will always be significant when x < 2, no matter how small 
z may be. 
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2.3. Solution at the critical point T = 0+ 

An important consequence of introducing this new form of interaction is that the system 
now becomes critical in the hmit T —>■ 0^ . This can be most clearly seen by considering 
a finite system of N subsystems with barriers 6j , i = 1 . . . N. For small T the system 
will remain static for long periods, but when a subsystem does eventually hop, the 
probability that it had the barrier bi is 



-b,/T 



Pi = ^-H.fT ■ (11) 



Since the system is finite it is always possible to identify a unique minimum barrier 
bmin (assuming that p{b) does not contain any delta-function peaks). Suppose it is 
subsystem j that has the barrier 6min , then inspection of (|Tl|) shows that pj — > 1 as 
T — s> 0+ while pt ^ for all i ^ j. Thus the identification of the activated subsystem is 
entirely deterministic - it is always the one with the smallest barrier. This is now 
algorithmically identical to the Bak-Sneppen model, already widely studied in the 
context of self-organised criticality |T^, |12|. Hence T = 0"^ is a critical point of the 
current model. Note that the timescale in the Bak-Sneppen model is normalised to 
precisely one hop per unit time, whereas in our model the time between hops can vary. 
However, as explained below, the maximum time between hops is of the order of e^''^'^ 
with be constant, which is finite for all T > 0. Hence we do not anticipate any problems 
with e.g. diverging moments of the distribution of hopping times. 

The Bak-Sneppen model is usually defined on a lattice, in which interactions only 
occur between nearest neighbours. This suggests that the current model can also be 
given a lattice interpretation, and this modification is considered in section ^. However, 
a mean field solution to the Bak-Sneppen model already exists [113, which should 



be the same as the T — 0^ solution of the current model. This indeed proves to be the 
case for the stationary distribution Poo{b), derived in [Appendix B| for arbitrary p{b), 



Pooib) ^ ^ ,/j!LyT ^--11 (12) 

~ ^(^ - ^c) Pib) as T ^ 0+, (13) 

where 9{x) is the usual theta function and be is defined by 

pip) db = ^ . (14) 
z + 1 

This is also the mean field solution to the Bak-Sneppen model, to within 0{1/N). 
Note that (|12D was found by taking T — > 0^ after taking the thermodynamic limit 
oo. That it agrees with the mean field solution of the Bak-Sneppen model, which 
corresponds to taking T — > O"*" before N ^ oo, suggests that the order of the limits 
may not be important for time-independent quantities, although this may not hold for 
time-dependent measures. 

The discontinuity in b in (|T3|) arises from the existence of a single characteristic 
time between interactions of the order of ujq e^"^^. The expected time until a subsystem 
with a barrier b hops is c^oe''/^, so subsystems with b <ti b^ will almost certainly hop 
before interacting, whereas those with b ^ b^ will almost certainly interact before 
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hopping. Only barriers b ^ be will have comparable probabilities of either hopping or 
interacting. As T is lowered the distribution of hopping times becomes extremely broad, 
corresponding to an increasingly narrow range of b for which the rates of hopping to 
interaction are comparable, eventually collapsing onto the single point be at T = 0+. 
Since hopping takes place over timescales much shorter than interactions, there will 
only be a vanishingly small proportion of subsystems with 6 < 6c at any given time, as 
implied by the theta function in (|1^). 

This suggests the classification of subsystems as either 'active' (those with b < b^) 
or 'inactive' (those with b > b^). In this heuristic scheme, active subsystems may become 
inactive by hopping into a trap with a barrier greater than be , but inactive subsystems 
can only become active via interactions. This description is somewhat reminiscent of 
the contact process, a critical model which falls into the universality class of directed 
percolation However, the contact process only becomes critical when its control 



parameter is set to some finite value - in effect, its value of be has to be set 'by hand.' 
By contrast, the current model automatically finds be in the limiting case T — O"*". It 
is in this sense that the Bak-Sneppen model can be said to be "self-organised" critical. 
A fuller discussion of the relationship between this model and self-organised criticality 
will be made elsewhere IJl^. Note that the Bak-Sneppen model and the contact process 



have different driving terms and fall into different universality classes ||. 
3. Finite dimensional lattice model 

The analysis of the preceeding section is mean field in the sense that the interactions 
have been assumed to act homogeneously throughout the system. In this section the 
model is given spatial definition by identifying each subsystem with a single site of 
a (i-dimensional lattice. The subsystems become thermally activated as before, but 
now each activation event can only alter the barriers of adjacent lattice sites, so the 
interactions are strictly short ranged. Only a one dimensional lattice is considered here, 
but the model can be defined in higher dimensions in a similar manner. 

The lattice model is defined as follows. A single barrier height bi is stored in each 
site i = 1 . . . N of a one dimensional lattice of size N, with initial conditions to be 
specified below. At the start of each time step, every site is checked to see whether it 
becomes activated, which it does with a probability . Any activated subsystems 

hop to a new trap with a new barrier bi drawn from p{b). Furthermore, the barriers 
in the adjacent sites 6i_i and 6j+i are also assigned new values - this represents the 
strong interactions. Note that in these scheme each site hops at most once per time 
step. This discrete time variable differs from the continuous time employed in the mean 
field equations in section |^, but should make little difference when the overall rate of 
activity is low. 

Thus every activation event causes three barriers to change value - the activated 
site itself, and the barriers in both of the adjacent sites. This would appear to fix 
z = 2. However, smaller values of z can be incorporated by stipulating that 6j_i and 
6j+i only change values with a probability z/2, where now Q < z <2. Simulation results 
for various z and T are discussed in sections (^TTj) to ( |3.3| ) below. Periodic boundary 



conditions have been assumed throughout, so b^^i = bi . The prior distribution was 

t The contact process also has shghtly different interaction terms to the Bak;-Sneppen model, but these 
can be altered without changing its universality class |l9| . 
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taken to be the simple exponential p{b) = ^e~^/''o. Care was also taken to ensure 
that the supplied random number generator gave good statistics deep in the tail of 
an exponential distribution. 



3.1. Existence of a stationary solution 

The simulations were first employed to see whether or not there exists a stationary 
solution in the lattice model, and if so, how robust it is with respect to the initial 
conditions. To test this, two very different initial configurations were used. The first 
was an initial 'quench' where all the bi were drawn from p{b) but spatially uncorrelated. 
These runs were then repeated with the same parameters but starting with every b^ 
arbitrarily large except for a single 'seed' at the origin with bo = (or any finite value). In 
both cases the barrier distribution P(6, t) appeared to converge to the same distribution 
Poo{b), which remained steady within the timeframe of the simulations, typically 10^ 
time steps (although the situation is not so clear for the T = 0"*" limit in low dimensions; 
see section \i.'S\) . Since the initial conditions were so very different, it seems likely that 
this same distribution is approached irrespective of the initial state. The stationarity of 
this solution is confirmed in section ^.2| . 

Although it has already been shown that the mean field equations of section 
admit a stationary solution, it is by no means clear that it is always, if ever, reached. To 
investigate this, the simulations were also repeated using a 'random nearest neighbour' 
(RNN) algorithm in which the neighbours of the hopping sites are chosen entirely at 
random from the remainder of the system, with new neighbours being chosen at every 
time step. Apart from finite size effects and the discrete time variable, this should 
behave identically to the mean field model. Simulations of the RNN model show that 
the same stationary solution was reached for both sets of initial conditions, as in the 
one dimensional case. Furthermore the numerical estimate of c^oo was found to agree 
well with the theoretical predictions (^ and (||), as plotted in Fig. |^. 

Also plotted in Fig. ^ are the equivalent values of c^oo for the one dimensional model, 
which consistently lie under the mean field values. This is because the interactions are 
now spatially correlated with the hopping sites and are more likely to occur in regions 
of high activity, ie. low barriers. Since the interactions tend to decrease barrier values 
towards the mean of p{b), their effect on P{b,t) will be weaker in the one dimensional 
model than the corresponding mean field system. Thus the stationary state will not be 
reached until much later times, when the activity will have decayed to smaller values, 
as in Fig. ^ In the T — 0^ limit the limiting activity is cjoo ~ (see [Appendix Bj ), 



so as a corollary b^ should be higher in lower dimensions, and indeed this has already 
been observed in the Bak-Sneppen model |T3 . 



3. 2. Interruption of aging 

An important feature of many glassy systems is that their behaviour depends upon the 
time since the initial quench, a property known as 'aging.' This phenomenom can be 
quantified by use of the two-time correlation function C(t„ + 1, tw)? defined here as the 
proportion of sites that have not hopped between times t^ and t^ + t following a quench 
at t = 0. In Bouchaud's model the system becomes time-translationally invariant above 
the glass transition, ie. C(t^ + t,t^) — >• Ceq(t) as t„ — oo when x > 1. By contrast. 
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Figure 2. Plot of Woo against x for z = 1 (upper data sets) and z = 0.1 (lower data 
sets). In both cases, the circles and squares correspond to simulation results for the 
random nearest neighbour and one dimensional models, respectively. The system size 
was N = 10** and the estimated error bars are smaller than the symbols. The solid 
lines are the mean field predictions from (R) and M). 



when X < 1 the system behaviour cannot be decoupled from t^ and instead one finds 
C(t„ + t,tw) — *■ Cneqit/tw), demonstrating that the system ages 0. 

To confirm that the simulations do indeed converge to a stationary solution, it 
is necessary to show that C(tw + ^,^w) converges to a time-translationally invariant 
solution Coo{t)- This has been checked in both one dimensional and RNN simulations, 
and in all cases stationarity was reached after waiting times t^ 2> tia , where depends 
on z, T and the dimensionality. A typical example is given in Fig. ^. An estimate 
of tia for low temperatures in the mean field can be found by following the procedure 
described in [0. This involves substituting the trial solution 



P{b,t) = ^u(f){u) , 

expressed here in terms of the dimensionless variable u 
the master equation (Rf) then becomes 



J_Qb/T 

LxJot 



For p{b) 



du 



'l/ujQt 




zu(p{u) — 



x{l + z)' 

[uq uty 



(15) 



(16) 



where the integral is just uj{t) expressed in terms of and u. Dimensional analysis 
suggests that the first term on the right hand side of (|16|) can be neglected for small t. 
This allows a t- independent scaling solution to be found which exhibits aging @]. 
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Figure 3. The two-time correlation function C(tw +t,tw) plotted against t for a one 
dimensional lattice with N = lO"*, z = 0.1, x = 0.35 and ujq = 1. From left to right, 
tw = 10", 10\ 10^, 10^, lO'^, 10^ and 10^. Note that the last two Hnes overlap. The 
system was started from an initial quench with the bi drawn from p{b) and uncorrelated. 

However, if 2; > then, as t — > 00, the second term will become negligible and it is 
straightforward to show that no self-consistent scaling solution exists. From this we 
infer that the aging is interrupted at some finite time. The crossover time tig, corresponds 
to times when both terms are of comparable magnitude, which can be estimated by 
dimensional analysis to be 

iia-(^)'^'~e^^/^ (17) 

where be is the same as that defined in (|1^. Note that tia diverges rapidly as T — 0, 
so aging behaviour may persist over intermediate (possibly experimental) timeframes at 
low temperatures. 

3.3. Aging in the T = 0+ limit 

The analysis of the T = Q"*" limit presented in section |2.3| was limited to the stationary 
barrier distribution Poo{b). We now turn to consider time-dependent behaviour, which 
requires careful consideration of the choice of timescale. Clearly as T — >^ O"*" the rate 
of hopping in a finite system becomes vanishingly small and so the timescale must be 
normalised in some manner to attain non-trivial behaviour. Two different timescales will 
be considered here. The first is the event timescale r in which precisely one subsystem 
from a total of hops per unit r. This is the usual choice in the Bak-Sneppen model. 
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The second time variable we shall consider is the avalanche timescale r^^, defined such 
that T^^ increases by one unit for every avalanche in the system, where an avalanche 
is defined below. Although not usually considered in connection with the Bak-Sneppen 
model, T^'^ is more in line with the timescale employed in e.g. the sandpile model pO[ . 



The existence of aging behaviour depends crucially upon which timescale is adopted, as 
we now demonstrate. 

Boettcher et al. have found that the Bak-Sneppen model exhibits a form of aging 
in one and two dimensional simulations starting from a single active seed We 
have adapted our code to employ event time in simulations starting from an initial 
quench and have also observed aging in the one dimensional model, but not in the RNN 
model. Results from the one dimensional simulations are given in Fig. |, which appears 
to exhibit aging for the longest times we were able to simulate, up to r„ = 10^ on a 
N = 2 X 10^ lattice, although the statistics become very noisy at later times. For large 
the data shows some indications of collapse onto a scaling function Coq{t/t^), as 
demonstrated in Fig. By contrast the results from the RNN simulations (not given) 
fail to exhibit any form of aging whatsoever. 

It is possible to provide a rough derivation of the aging behaviour observed in our 
simulations by considering the pattern of activity of hopping sites as the system evolves 
from its initial quench. This argument employs known results for the Bak-Sneppen 
model, which we do not attempt to justify here; the reader is instead referred to 



for a complete description. Initially the location of the hopping site moves around 
the system at random, but as r increases it tends to repeatedly visit one part of the 
system for a finite time before jumping to another, uncorrelated region. These localised 
spatiotemporal regions of activity are known as avalanches. The expected duration of 
the avalanche (S) and the expected number of sites covered during the avalanche (ucov) 
diverge with time according to |T2| 



1 

r \ 7-1 



Kov) ~ [-^) , (19) 

where the critical exponent 7 depends on the dimensionality. In low dimensions 7 > 1, 
but 7 = 1 in high dimensions where the scaling relations given by (|THp and (|TU|) break 
down. In this case the active site never remains localised in any one part of the system 
for any significant time and the following argument does not apply. 

As before, let C{t^ + t,t^) be the proportion of the system that has not hopped 
between times and + r. Clearly C will decrease when the avalanche jumps to a 
new part of the system. If we now coarse grain over length and time scales much larger 
than {S) and (ucov), then, as consecutive avalanches are spatially uncorrelated, we can 
write the following continuous differential equation for C, 

dC _ C (ncov) 
dr^ N {S) • 

C only decays if the avalanche jumps to a part of the system it has not already visited. 



hence the factor of C on the right hand side of (pOD . Furthermore each avalanche changes 
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a proportion {ncov)/N of the system and lasts for a time (S), which accounts for the 
remaining factors. Integrating (pO) from to + r with C{t^, t^) = 1 gives 



(21) 



C(rw + r, Tw) ~ (^1 + 



where a = \imr-,oo(yT {ucov))/ {N {S)). That C(r„ + r, r^) depends on r and only 
through the ratio r/r^ demonstrates that the system ages. Although this derivation 
is sufficient to explain our simulation results, it is not clear that the coarse graining 
assumption is valid as r ^ oo, even for arbitrarily large systems. It is also not obvious 
if this procedure can be extended to account for the results of Boettcher et al. |22 



If indeed the Bak-Sneppen model does age in event time r, this calls into question 
the usual assumption of stationarity in most previous studies of the Bak-Sneppen model, 
since an aging system does not obey time-translational invariance and hence cannot be 
stationary. It is possible that previous studies have relied too heavily on one-time 
functions such as critical exponents etc., which may appear to tend to stationary values 
even when two-time functions are still evolving. This is a subtle question and further 
investigation would be desirable. Note that the aging we have observed (apparent or 
otherwise) is dimensional-dependent — it only holds in low dimensions, not in the mean 
field. Hence the assumption of stationarity in section |2.3| is still valid. 

Turning now to consider avalanche time r^^, the same coarse-graining assumptions 
lead to the following differential equation for C{t^ + t^^ , r^^). 



dC 

which differs from ( 



N ' 

by a factor of (5). This can be solved as before to give 



,rr)~exp<!-a(^ 



7-1 



(22) 



(23) 



where a is an arbitrary constant. Since this expression tends to zero as oo for all 

> 0, we conclude that the system does not age in avalanche time. This may relate to 
the failure to observe aging in the sandpile model |2l| , where a timescale similar to r^^ 



is usually employed. Note that this conclusion also holds in the mean field case 7 = 1, 
although the form of (|23|) is different. The numerical plot of C{t^ + r^^, t^) for a one 
dimensional system is given in Fig. 0. 



4. Conclusions 

We have introduced the strongly-interacting trap model, a version of Bouchaud's 
trap model which exhibits glassy behaviour for low temperatures and intermediate 
timeframes, but ultimately reaches stationarity as t ^ oo. In the limit of zero 
temperature the model becomes critical in all dimensions, including the mean field, 
so it has no 'lower critical dimension.' 

A lattice version of the model was also introduced which appears to behave 
qualitatively similar to the mean field model in all the cases we have looked at, except 
one. This is the apparent aging in the zero-temperature limit in low dimensions in what 
we have termed 'event time,' which is absent in the mean field. This would represent 
a form of dimensional-dependent aging, but is currently only supported by numerical 
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Event time x 



Figure 4. C(tw + t, Tw) versus r for a one dimensional lattice of -/V = 2 x 10''' sites 
in the T — > 0'*' limit in event time r, making this identical to the Bak-Sneppen model. 
From left to right, the different lines refer to t„ ^ 10°, 1Q\ 10^, 10^, 10"*, 10^ 10^, lO'^ 
and 10^, respectively, although the first 4 lines overlap. The apparent cut-off for small 
Tw is purely the result of an optimisation procedure employed to improve run times. 

evidence and hence remains somewhat speculative. It is also not clear at what dimension 
the aging behaviour might start to break down. 

Although this work was originally motivated by the range of spin glass models that 
have a critical point at T = 0, we do not claim any precise relationship between these 
systems and our model. We merely propose that our model may serve as a caricature 
of glass models with zero-temperature criticality, which, by its very simplicity, should 
allow for fuller analysis of this class of systems, both theoretically and numerically. It 
may also serve as a link between glass theory and models of self-organised criticality. 
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Appendix A. Derivation of (0) and (§) 

In this appendix we derive the expressions for quoted in (|^) and (H) for non-integer 
and integer values of x = T/bo , respectively. The steady state solution is found by 
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O 
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0.4 




10' 

Event time x/t„ 



10^ 



Figure 5. The same data as in Fig. || for Tw > 10^, this time plotted against t/tw , 
indicating partial collapse. From right to left (as viewed from the upper portion of the 
graph), the lines refer to Tw — 10^, 10^, 10'' and 10^ respectively. 



substituting the form of Poo{b) given by (|^) into the definition of u{t) given in (Q). This 
results in the following integral equation for , 

1 =f^^db. (A.1) 

Substituting p(6) = ^e^''^^" and making the change of variables u = e'^^^ gives 

X / du (A.2) 



z+1 Jo U + ZUJ, 

= ^—-F(i,1+x-2 + x;-{zuj^)-^). (A.3) 

For z ujoo < 1 and non-integer x the hypergeometric function F can be rewritten as [B^ 



,ZUJr 



ZUr, 



X/ T(2 + x)T(x) ^ i-zuj^ 

T{2 + x)T{-x){zUooT -X /JJ Y. 



r(l + a;))2 



k=0 



k — X 



(A.4) 



where T{x) is the usual gamma function. Inserting this into ( [A. 3D gives the final 
expression (0), where use has been made of the identity 



Tix]T{-x) 



X sm Tlx 



(A.5) 
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Figure 6. C{t^ + t'^^jT^) measured in avalanche time t'^^ for a one dimensional 
lattice with N — 10"* sites in the T = 0+ limit. From upper-right to lower-left, the lines 
refer to r^^ = R^, R' , R^, R^ and R^° respectively, where R = 2.512. 



The sum on the right hand side of ( |A.4| ) is ill defined when x is an integer. In this 



case, an equivalent expression can be found, either by integrating ( [A. 21 ) directly or by 



substituting x = n^e into (|^) with n an integer, and letting e Q. Both methods 
yield the same answer, given in (^). 

Appendix B. Derivation of the T = 0+ solution ( [T^ 

It is possible to derive, in a non-rigorous manner, explicit expressions for u^o and Pooip) 
for arbitrary p{h) in the limit T O"*". The starting point is the same integral equation 
as in ( |OD , 

1 - db . (B.l) 



Z + \ Jo Qb/T+ln{zu}oa) _|_ X 

This imposes significant constraints on the allowed forms of Uoo ■ For instance, if 
T\ia{z ujoo) — >■ as T — >• 0, then (p.l| ) would imply the contradictory result z = oo. 
Indeed, if one assumes that \ia{zu!oo) ~ T"", then self-consistency of (p3.1| ) demands that 
a = — 1 and therefore 

^oo ~ - e-"^/^ , (B.2) 



where be is some arbitrary constant. 
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For \n{zujoo) = —hc/T, the denominator inside the integral of ( [B.l| ) is essentially 
constant except around the region b ^ be ■ Since p{b) is independent of T, we can always 
choose T small enough so that p{b) is slowly varying over this region, assuming that 
p{b) is continuous at be ■ Hence the integral simplifies to 

= p{b) db , (B.3) 

z + I Jo 

which fixes be ■ Substituting ( |B.2|) into the general expression for Poo(^) given in (|^) 
results in the final expression (|1^). To give some idea of typical values of 6c , a simple 
exponential prior p(6) = ^e~^/''o gives be = bQ\n{l + z~^), whereas be = {1 + z)^^ when 
pip) is uniformly distributed on [0, 1]. 
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